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SECTION  I 
INTRODUCTION 

This  report  is  concerned  with  the  refinement  of  the  computer  code 
SINGER,  which  was  developed  under  contract  F29601 -73-C-0022.  The  function 
of  SINGER  is  to  predict  the  nonlinear  response,  including  member  failures 
and  structural  collapse,  of  plane  skeletal  reinforced  concrete  structures 
to  static  and  dynamic  loads.  The  nonlinearities  may  be  caused  by  finite 
displacements  and  inelastic  constitutive  laws  [ 1]. 

The  mathematical  model  of  SINGER  represents  an  assemblage  of  one- 
dimensional finite  elements,  which  simulate  the  behavior  of  beam-columns. 

The  model  is  expressed  in  the  form  of  an  energy  function  which  is  defined 
by  the  generalized  coordinates  and  forces  (external  and  inertia).  The 
origin  of  the  generalized  coordinates  corresponds  to  the  unstrained  state 
of  the  system,  the  initial  state. 

The  motion  (equilibrium  path)  of  the  system  is  determined  at  a dis- 
crete number  of  points  in  time  (load  increments).  The  solution  is  ex- 
tended through  successive  time  steps  by  the  repeated  application  of  the 
solution  process,  which  comprises  the  following  tasks: 

1.  The  state  of  the  system  is  defined  at  the  beginning  of  a time  step. 

2.  The  motion  during  the  time  step  is  discretized  such  that  the  ac- 
celeration at  the  end  of  the  time  step  is  defined  in  terms  of  the 
initial  conditions  and  the  unknown  generalized  coordinates  at  the 
end  of  the  time  step. 

3.  The  energy  funtion  is  formulated  at  the  end  of  the  time  step  and 
minimized  to  yield  the  generalized  coordinates. 
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It  follows  that  the  solution  contains  spatial  and  temporal  discreti- 
zation errors.  Temporal  error  controls  are  contained  in  SINGER  [ 1]. 
Spatial  error  controls  are  proposed  in  this  report. 

PURPOSE  AND  SCOPE 

The  refinement  of  SINGER  is  based  on  the  following  principal  tasks 
(sections  II  - IV):  The  validation  of  the  finite  element,  the  selection 

of  an  appropriate  quadrature  method,  and  the  formulation  of  error  controls 
in  energy  evaluations.  These  studies  motivated  the  development  of  a new 
finite  element,  which  is  presented  in  section  V.  The  current  capability 
of  SINGER  is  demonstrated  in  section  VI,  and  modifications  of  SINGER  (im- 
plemented and  proposed)  are  discussed  in  section  VII. 
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SECTION  II 


VALIDATION  OF  THE  FINITE  ELEMENT  MODEL 

This  section  is  concerned  with  the  validation  of  the  macro  response  of 
the  finite  element  model  of  a beam-column. 

FINITE  ELEMENT  MODEL 

Beam-column  models  represent  slender  elements  whose  axial  and  flexural 
deformations  are  geometrically  coupled.  The  finite  element  model  of  a beam- 
column,  the  basic  component  of  the  discrete  model  in  SINGER,  is  based  on 
the  following  assumptions  of  the  classical  beam-column  [2,  3J: 

1.  The  plane  of  bending  coincides  with  the  longitudinal  plane  of 
symmetry. 

2.  Plane  sections  remain  plane  and  normal  to  the  deformed  reference 
axis. 

3.  Normal  strains  and  rotations  are  infinitesimally  small. 

4.  The  stress-strain  law  of  any  lonqitudinal  fiber  is  defined  by  the 
constitutive  law  of  the  material. 

5.  The  effect  of  shear  deformations  is  negligible. 

These  assumptions  have  been  substantiated  for  elastic  and  inelastic  reinforced 
concrete  and  steel  beam-columns,  e.g.,  [ 2,  4 , 5].  Accordingly,  the  classical 
beam-column  provides  the  basis  for  validating  the  finite  element  model. 

It  should  be  noted  that  in  spite  of  the  restriction  of  the  finite  ele- 
ment model  to  small  deformations,  there  are  no  limitations  on  the  displace- 
ments of  the  nodes  to  which  the  elements  are  connected  [ g ] . 

VALIDATION  BASIS 

The  effectiveness  of  the  finite  element  model  is  judged  on  the  basis  of 
its  performance  relative  to  classical  continuum  models  of  columns  and  beam- 
columns.  The  performance  of  the  model  is  measured  in  terms  of  micro  and 


* Numbers  in  brackets  designate  references. 
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macro  parameters.  Micro  parameters  consist  of  strains,  stresses,  and 
internal -energy  densities;  macro  parameters  describe  loads,  displace- 
ments and  element  energies. 

This  distinction  in  model  parameters  is  made  to  reflect  the  difference 
between  pointwise  and  average  behavior  of  discrete  models.  For  instance, 
Oliveira  [ 6]  proved  that  the  potential  energy  of  an  assembly  of  complete 
and  conforming  finite  elements  converges  in  the  limit  to  the  potential 
energy  of  the  corresponding  continuum  model.  However,  this  does  not  guaran- 
tee uniform  (pointwise)  convergence  of  strains  and  stresses.  Similarly,  in 
the  Ritz  method,  the  convergence  of  a sequence  of  approximating  functions 
to  the  exact  solution  can  only  be  guaranteed  in  a mean-square  (average) 
sense  [ 7]. 

VALIDATION  PROBLEMS 

The  performance  of  the  finite  element  model  (and  assemblages  of  the 
model)  is  tested  over  the  entire  range  of  response,  from  the  elastic  geome- 
trically linear  range  to  the  inelastic  geometrically  nonlinear  range. 

In  this  section,  the  macro  response  predictions  of  the  discrete 
element  (and  system)  model  are  compared  with  those  of  two  continuum  models: 
the  Elastica  model  [ 8],  and  the  Jezek  model  [9].  These  continuum  models 
have  been  selected  for  the  validation  purpose  because  they  admit  closed- 
form  solutions.  The  evaluation  of  the  micro  response  is  conducted  in 
sections  IV  and  VI. 

The  Elastica  problem  was  first  investigated  by  Euler.  (A  historical 
discussion  is  presented  by  Timoshenko  [10]).  It  is  concerned  with  the  so- 
lution of  the  'exact'  differential  equation  of  a buckled  column  for  the 
elastic  curve  (the  deformed  reference  axis).  The  term  exact  refers  to  the 
formulation  of  the  curvature  of  the  elastic  curve.  Accordingly,  assumption  3 
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of  the  classical  beam-column  model  is  relaxed.  However,  the  remaining 
assumptions  are  inherent  in  the  Elastica  model.  The  closed-form  solution 
of  the  Elastica  model  is  obtained  in  the  form  of  elliptic  integrals  [ 8 J - 
The  Jezek  model  [9]  was  developed  to  obtain  a rigorous  analytical 
solution  to  the  stability  problem  of  inelastic  beam-columns.  Karman  [3] 
was  the  first  [ 2 ] to  investigate  the  buckling  of  an  inelastic  eccentrically 
loaded  column  as  a stability  problem.  Karman  formulated  the  differential 
eguation  governing  the  eguilibrium  path  of  the  column  consistent  with  as- 
sumptions 1 through  5 of  the  classical  beam-column.  In  addition,  he  as- 
sumed that  no  strain  reversals  occurred  during  any  stage  of  the  loading 
process,  including  the  attainment  of  the  limit  point.  The  constitutive 
law  was  obtained  from  a unidimensional  compression  test.  (The  eccentri- 
cities were  small  enough  to  preclude  to  apprarence  of  tensile  strains). 

The  solution  was  obtained  with  great  precision  by  numerical  integration. 
Karman's  theory  and  his  results  were  validated  by  Chwalla  in  a series  of 
theoretical  and  experimental  investigations  [5,11] 

Jezek  [ 9]  simplified  the  Karman  model  by  introducing  an  elastic- 
plastic  stress-strain  law.  This  permitted  him  to  divide  the  beam-column 
■>nto  subregions  of  at  most  three  distinct  stress  states;  I.  elastic,  II. 
inelastic  in  compression,  III.  inelastic  in  compression  and  tension  . He 
formulated  differential  equations  for  the  distinct  stress  states  (2nd  order 
differential  equations  in  the  transverse  deflection,  which  are  nonlinear 
for  states  II  and  III)  and  obtained  closed-form  implicit  relations  between 
the  load  parameter  and  the  transverse  deflection.  The  constants  of  in- 
tegration were  evaluated  by  imposinq  boundary  conditions  and  continuity 
conditions  at  the  interfaces  of  the  subregions.  He  presented  solutions 
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for  a rectangular  beam-column  with  two  loading  conditions;(l ) eccentri- 
cally loaded  and  (2)  concentrically  loaded  with  a concentrated  transverse 
load  at  the  midspan. 

Recently, closed-form  solutions  to  Jezek's  problems  were  presented 
again  [12,13],  They  differ  from  Jezek's  work  as  follows:  The  curvature 

rather  than  the  transverse  deflection  is  selected  as  the  dependent  vari- 
able, complex  as  well  as  rectangular  cross  sections  are  considered,  and 
the  axial  loads  of  the  beam-column  with  the  concentrated  load  at  the  mid- 
soan  can  be  applied  eccentrically.  Aside  from  that,  the  models  and  the 
solution  procedures  of  the  two  versions  are  identical.  These  publications 
resulted  in  a computer  program  [14]  which  is  used  in  the  validation  effort 
of  inelastic  beam-columns. 

ELASTIC  RESPONSE 

A detailed  study  of  the  macro  response  of  the  discrete  model  in  the 
elastic  geometrically  linear  and  nonlinear  range,  with  special  attention 
to  the  performance  of  the  internal  node,  is  presented  in  reference  15. 

A summary  of  the  major  results  of  this  successfully  completed  validation 
phase  is  presented  below. 

The  function  of  the  internal  node  is  to  make  the  location  of  the  re- 
ference axis  independent  of  the  strain  state  [1  ].  It  was  decided  [1  ] that 
this  goal  can  be  accomplished  if  the  finite  element  is  capable  of  represent- 
ing a linearly  varying  normal  strain  along  its  reference  axis.  For  this 
purpose  the  internal  node  was  introduced. 

This  study  indicates  that  the  internal  node  does  successfully  perform 
its  function.  The  response  predictions  (micro  and  macro)  are  independent 
of  the  location  of  the  reference  axis;  i.e.,  the  reference  axis  may  be 
chosen  to  coincide  with  any  longitudinal  fiber  of  the  beam-column  or  may 
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even  be  placed  outside  the  beam  without  altering  the  response  predictions. 

A postbuckling  analysis  of  a concentrically  loaded  elastic  column  is 
conducted  to  evaluate  the  ability  of  the  finite  element  model  to  predict 
large  displacements  with  high  precision.  The  results  of  the  analysis  are 
presented  in  figures  1 and  2 for  assemblages  consisting  of  4 and  9 elements, 
respectively.  They  indicate  that  a sequence  of  finite  element  approxima- 
tions converges  monotonical ly  to  the  exact  solution  of  the  Elastica  model 
and  that  any  desired  degree  of  accuracy  can  be  attained  by  refining  the 
element  mesh. 

INELASTIC  RESPONSE 

The  second  Jezek  problem  [9],  a concentrically  loaded, simply  supported 
beam-column  with  a concentrated  transverse  midspan  load  (figure  3a),  forms  the 
basis  of  the  validation  effort  in  the  inelastic  geometrically  nonlinear  range. 
The  solution  of  the  continuum  model  is  obtained  via  the  Lehigh  program  [14] 
for  various  load  combinations  (for  each  prescribed  axial  load  the  transverse 
load  is  varied  until  the  limit  load  is  attained)  and  finite  element  mesh 
sizes.  The  response  is  represented  by  the  equilibrium  path  ( the  midspan 
load  vs.  midspan  deflection)  and  compared  with  that  of  the  discrete  model. 

The  results  presented  in  figure  4 are  characteristic  of  the  test  pro- 
blems studied:  The  finite  element  solution  expressed  by  the  equilibrium 

path  converges  to  the  exact  solution,  and  the  accuracy  of  the  macro  response 
can  be  controlled  via  the  mesh  size. 

In  the  course  of  this  study  it  became  apparent  that  the  finite  element 
model  has  a tendency  to  develop  a directional  preference  with  increasing 
deformations.  For  instance,  when  the  1-axes  of  all  elements  point  from 
left  to  right  (fiqure  3a),  a slight  discrepancy  in  the  transverse  deflect 
ions  of  nodes  2 and  4 appears  at  small  loads  and  grows  with  the  load  leve 
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FIG.  4 LEHIGH  COMPARISON  PROBLEM 
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On  the  other  hand,  when  the  elements  are  symmetrically  arranged  (the  1- 
axes  of  all  elements  point  towards  the  center  as  in  figure  3b),  the  trans- 
verse deflections  of  nodes  2 and  4 are  identical. 

The  directional  tendency  of  the  finite  element  model  appears  to  result 
from  a combination  of  discretization  and  roundoff  errors  and  is  magnified 
by  the  local  frame  of  reference,  which  does  not  reflect  the  inherent  sym- 
metry of  the  element  geometry.  The  discretization  error  is  associated  with 
the  approximate  representation  of  the  deformation  field  (and,  hence,  the 
strain  and  internal -energy  density  field)  of  the  continuum  model  by  four 
interpolation  functions  corresponding  to  the  deformation  components  u^u^, 
0^,  (T4  (rigid-body  motions  have  been  eliminated).  The  components  lT-|  , u,?,  u3 
are  obtained  via  coordinate  transformations  from  the  global  displacements 
of  the  nodes  to  which  the  element  is  connected.  Figure  3c  indicates  that 
for  the  same  nodal  displacements,  one  can  obtain  two  distinct  sets  of  de- 
formation components  depending  on  the  orientation  of  the  local  1-axis 
(left  to  right,  or  right  to  left).  Since  each  set  contains  inevitable 
roundoff  errors,  the  corresponding  deformation  fields  will  also  contain 
roundoff  errors  and  are  not  likely  to  be  identical.  This  accounts  for  the 
directional  preference  of  the  element. 

To  minimize  the  directional  tendency  of  the  finite  element,  a new 
model  has  been  developed  based  on  a local  frame  of  reference  that  is 
intrinsic  to  the  element  geometry  (see  section  V).  This  model  contains 
also  a refinement  in  the  representation  of  the  strain  field,  which  was 
motivated  by  the  study  of  error  controls  (section  IV). 
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SECTION  III 


QUADRATURE  METHOD* 

This  section  is  concerned  with  the  selection  of  an  integration  method 
to  compute  the  internal  energy  of  the  finite  element  with  sufficient  ac- 
curacy to  assure  energy  convergence.  The  internal  energy  of  the  element 
is  defined  by  the  relation 

n = J IT*  dV  (1 ) 

V 

where  the  internal -energy  density 
* e 

n = f ade  (2) 

0 

In  Eqs.  1 and  2,  V represents  the  volume  of  the  element,  and  a and  e 
denote  the  stress  and  strain  at  any  point  in  the  volume,  respectively. 
Energy  evaluations  must  be  performed  numerically  because  of  the  complex 
variation  of  the  internal -energy  density  over  the  volume  of  a nonlinear 
beam-column.  A reliable  method  that  requires  the  least  effort  to  produce 
a specific  degree  of  accuracy  (i.e.,;  the  "best"  quadrature  method)  is 
selected  to  evaluate  Eq.  1.  A detailed  literature  search  shows  that  the 
Gaussian  quadrature  method  is  the  "best"  quadrature  method;  the  selection 
is  based  on  convergence  characteristics  and  error  estimates  of  several 
quadrature  methods. 

There  are  two  basic  procedures  for  evaluating  numerically  the  integral 
of  a function  of  one  variable  [16,17].  In  the  first,  points  at  which  the 


* Prepared  by  Ayodele  Abatan,  graduate  student,  VPISSU. 
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function  is  evaluated  are  specified  a priori  at  equal  intervals,  and  a 
polynomial  is  passed  through  these  function  values  and  integrated  exactly. 
In  the  second,  the  sampling  points  are  located  so  as  to  achieve  maximum 
accuracy.  These  two  procedures  lead  to  a broad  classification  of  quadra- 
ture methods  into  the  following  categories: 

1 EQUIDISTANT  interpolation-point  methods 
2.  VARIABLE  interpolation-point  methods 

EQUIDISTANT  INTERPOLATION-POINT  METHODS 


Quadratures  possessing  equidistant  interpolation  points  are  characterized 


by  the  Newton-Cotes  quadrature.  The  quadrature  formula  can  be  expressed  as 

n 

I = [fU)  = l H.  f(ef)  + E (3) 

i = l 


where 


£.j  = interpolation  points  equidistantly  spaced  over  the  normalized 
domain  -1  < c < 1 
= weights 

E = error  of  order  (a)11 


n = number  of  interpolation  points 

Eq.  3 is  obtained  by  approximating  the  function  f ( ^)  with  a Lagrange  interpo- 
lation formula  and  evaluating  the  resulting  integral  exactly.  For  n = 3, 

Eq.  3 yields 

I = J [f(-l)  + 4f (0 ) + f(l )]  (4) 

which  is  known  as  Simpson’s  rule.  Modifications  of  Simpson's  rule  and 
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other  rules  of  equidistant  interpolation-point  methods  are  presented  in 
the  literature,  e.g.,  [18]. 

VARIABLE  INTERPOLATION-POINT  METHODS 

The  variable  interpolation-point  methods  are  known  as  Gaussian 
quadratures.  Gauss  had  the  ingenious  idea  [18]  to  let  the  interpolation 
points  be  variables.  Hence  for  n interpolation  points  (Gauss  points), 
there  are  2n  unknowns  [c^ , f(s^),  i = 1,  2,  ...  n];  which  means  that 
an  n-point  Gaussian  quadrature  formula  is  exact  for  a polynomial  of 
degree  2n-l . Gauss  found  that  by  selecting  Legendre  polynomials  as 
interpolation  formulas,  he  could  solve  the  problem  in  explicit  form. 
Gaussian  quadrature  formulas  have  also  been  developed  for  other  inter- 
polation formulas  [19]  such  as  Hermite,  Laguerre,  and  Jacobi  interpolation 
formulas.  The  Gaussian-Legendre  quadrature  formula,  a common  and 
highly  convergent  formula,  can  be  expresstd  in  the  form 
1 n 

I = f fU)d£  = l H.  f(?i)  + E (5) 

-1  1=1 

where 

£.  = Gauss  Points 
H.j  = weights 
E = error 

The  location  of  Gauss  points  together  with  the  associated  weights 
are  tabulated  in  the  literature,  e.g.,  [19,  20,  21]. 
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The  quadratures  can  also  be  formulated  for  multiple  integrals: 

11  n m 

I = J J fU,n)d£dn  = £ I H.Hj.f (5.  ,nj ) + E (6) 

-1  -1  1=1  j=l 

where  n,  m represent  the  number  of  Gauss  points  in  the  £,n  direction, 
respectively. 

COMPARISON  OF  QUADRATURE  METHODS 
The  comparison  is  based  on 

a.  convergence  characteristics 

b.  error  estimates 

Equidistant  interpolation  does  not  in  general  lead  to  a well -convergent 
process.  Convergence  is  not  even  guaranteed  for  functions  which  have  singu- 
larities inside  the  unit  circle  in  a complex  plane  even  though  they  may  be 
smooth  in  such  a normalized  domain.  Convergence  is  more  likely  for  func- 
tions which  are  analytic  (i.e.,  the  functions  and  their  derivatives  are 
continuous)  within  the  solution  domain  [18], 

The  convergence  of  the  Gaussian  quadrature  process  is  guaranteed  by 
the  general  nature  of  orthogonal  expansions  [18].  Many  convergence  studies 
and  error  estimates  of  the  Gaussian  quadrature  are  based  on  the  definite 
integral 

b 

j w(x)  f(x)dx  (7) 

a 

where  w(x)  is  a function  for  which  the  moments  (or  nominal  integrals) 
b k 

Ck  = J w(*)xax,  k = 1 , 2,  . . . (8) 

a 

are  defined  and  C > 0. 
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Gaussian  formulas  of  the  form  qiven  by  Eq.  5 for  nonnegative  weight 
functions  have  several  important  convergence  properties.  These  properties 


are  the  fundamental  reasons  why  Gaussian  formulas  are  so  important  in 
applications.  These  formulas  converqe  to  the  true  value  of  the  integral 
for  almost  any  conceivable  function  which  one  meets  in  practice.  Also 
the  convergence  properties  for  analytic  functions  show  that  essentially 
the  only  formulas  which  are  interpolatory  (i.e.  exact  for  all  polynomials 
of  degree  <_  n-1  [20])and  which  converqe  for  all  functions  analytic  in  a 
region  of  the  complex  plane  containing  [a,  b]  are  the  Gaussian  formulas. 

Ivanova  [22]  presents  results  concerninq  the  convergence  of  Gaussian 
formulas  on  an  infinite  segment.  Krylov  [20]  discusses  convergence  of 
general  quadrature  formulas  on  a finite  segment.  The  following  theorem 
due  to  Stroud  and  Secrest  [19]  is  stated  here  without  proof: 

Let  w(x)  be  nonnegative  on  [a,  b]  and  let  f(x)  be  continuous 

on  this  segment.  Let  x^^and  H.^n\  i = 1,  2 n,  be 

points  and  coefficients  of  the  n-point  formula  of  deqree  2n -1 
for  w(x)  on  [a , b] , then 


n b 

lim  l H.n  f ( x . = f w(x)  dx 

i=l  1 1 \ 


Gaussian  formulas  converge,  however,  for  a wider  class  of  functions  than  the 

li 

class  of  continuous  functions.  Proofs  of  this  statement  are  qiven  by  Szego 
[23]-  Several  definite  statements  have  been  made  reqardinq  the  usefulness 

of  the  Gaussian  quadrature  method  [ 17,18,21] .For  example,  Lanczos  [18] 

states:  "In  engineering  applications,  it  happens  rather  frequently  that  the 

average  value  of  a function  of  unknown  structure  has  to  be  established  on  the 

basis  of  very  few  observations.  In  this  case,  it  is  strongly  advocated  that 
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the  points  where  the  ordinates  are  measured  shall  follow  the  Gaussian 
pattern.  For  example,  the  pressure  tubes  in  an  airduct  will  give  much 
more  favorable  results  if  they  are  not  uniformly  distributed  over  the 
cross-section  of  the  airduct  but  in  conformity  with  the  Gaussian  zeros. 
The  same  holds  for  temperature  measurements  along  a wall  if  the  purpose 
of  these  measurements  is  to  establish  average  values." 

The  error,  or  remainder,  E(f),  of  the  quadrature  formula  is  given  by 


b n 


Usually  error  estimates  [19]  are  provided  in  the  form 

|E(f)|  < C M(f)  (11) 

where  C is  a constant  depending  only  on  the  quadrature  formula  and  M(f) 
is  a bound  on  a quantity  related  to  f(x).  For  this  discussion,  M(f)  is 
regarded  as  a bound  on  the  rth  derivative  of  f(x).  Replacing  the  con- 
stant C by  e , Eq.  11  assumes  the  form 

| E(f ) I < er  M(f)  (12) 

The  constants  er  which  depend  on  the  quadrature  furmulas  and  which  must  be 
available  to  use  the  estimates  are  given  by  Stroud  and  Secrest  [l Q ] . For 

example,  for  a Gaussian  formula 
b 

er  = l<  J w(x)  [Pn(x)]2  dx,  1 <r  <2n  (13) 

a 

where  Pn(x)  is  a polynomial  with  the  leading  coefficient  equal  to  unity, 
whose  zeros  are  the  interpolation  points  of  the  quadrature  formula. 


22 


It  has  been  shown  on  the  basis  of  Eq.  13  that  not  only  are  the 
Gaussian  formulas  "best"  for  integrands  of  higher  order  derivatives,  but 
they  are  also  quite  accurate  for  functions  having  only  low  order  deriv- 
atives. Stroud  and  Secrest  [Id]  showed  that  for  r = 1,  2 the  Gaussian 
formulas  are  as  good  as  the  best  formulas  which  have  an  equal  number  of 
points . 

On  the  basis  of  the  comparison  of  the  Gaussian  quadrature  method 
with  equidistant  methods,  it  is  evident  that  the  Gaussian  quadrature 
method  excells  not  only  in  regard  to  its  convergence  characteri sties  but 
also  with  respect  to  the  associated  error  estimate  for  a specified  accuracy. 
Hence, the  Gaussian  quadrature  method  is  selected  as  the  "best"  quadrature 
method . 
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SECTION  IV 

ERROR  CONTROLS  IN  ENERGY  EVALUATIONS 

The  accuracy  of  the  response  predictions  via  SINGER  is  limited  by 
the  accuracy  inherent  in  energy  evaluations.  The  search  process  for 
equilibrium  states  is  governed  by  variations  in  energy  resulting  from 
changes  in  the  generalized  coordinates  of  the  finite  element  model. 

If  these  energy  variations  do  not  accurately  reflect  the  state  of  the 
system,  the  response  predictions  are  likely  to  be  of  poor  quality  and 
the  solution  process  may  even  fail  to  converge. 

There  are  two  types  of  errors  that  must  be  controlled  to  assure  accurate 
energy  evaluations;  they  are  discretization  errors  and  quadrature  errors. 

The  discretization  error  is  caused  by  the  representation  of  the  continuum 
model  of  a skeletal  reinforced  concrete  structure  by  an  assembly  of  finite 
elements.  The  quadrature  error  is  caused  by  the  numerical  integration  of 
the  internal -energy  density  over  the  volume  of  the  finite  element. 

In  order  to  control  the  discretization  error,  it  is  necessary  to  know 
the  convergence  characteristics  of  a sequence  of  finite  element  approxi- 
mation. The  finite  element  model  of  the  beam-column  satisfies  the  condi- 
tions of  completeness  and  conformity,  which  guarantees  monotone  converqence 
(in  energy)  for  linear  elastic  systems  [ 6j.  Consequently,  in  the  linear 
range  of  response,  any  refinement  of  the  finite  element  mesh  improves  the 
quality  of  the  energy  prediction  of  the  discrete  model.  However,  no 
general  convergence  criteria  are  available  for  the  inelastic  range.  In 
this  region,  converqence  characteristics  can  only  be  obtained  by  numerical 
experimentation . 

The  quadrature  error  can  be  controlled  directly  by  increasing  the 
number  of  integration  points  and  indirectly  by  refining  the  finite  element 
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mesh.  If  the  constitutive  law  is  expressed  in  polynomial  form,  the 
integration  points  can  be  selected  to  achieve  exact  quadratures. 
Otherwise,  'exact'  integration  point  distributions  can  only  be  obtained 
by  numerical  experimentation. 


DISCRETIZATION  ERROR 

A test  problem  has  been  selected  to  study  the  convergence  character- 
istics of  the  finite  element  model  over  the  entire  range  of  response  and 
to  establish  guidelines  for  discretization  error  control . The  test  problem 
consists  of  a single  element  fixed  at  one  end  and  subjected  to  a concen- 
trated moment  at  the  free  end.  The  percent  deviation  of  the  extreme  fiber 
strain  in  compression  at  the  fixed  end  from  the  exact  value  is  selected  as 
a measure  of  the  discretization  error;  i.e., 


ed 


c - c 100 


(14) 


where  ed  is  the  discretization  error  measure,  ec  is  the  exact  value, and  ec 

is  the  approximate  value  of  the  maximum  compressive  strain  at  the  fixed  end. 
From  Eqs.  2.12  to  2.19  in  reference  1 , ?c  can  be  expressed  as 


- _ a , b 

e„  - £„  + e_ 

c c c 


(15) 


where 

c*c  = (4u4  - iTj  )/L  (16) 

and 

ec  = "(6iVL  ‘ 2u^)/L  (17) 

The  reference  axis  coincides  with  the  centroidal  axis  of  the  cross-section, 
which  is  a 2 inch  square. 
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In  Eqs . 16  and  17,  lT,  , represent  axial  displacements  and  u2 , 
represent  bendinq  displacements. 

This  test  problem  has  been  selected  from  a group  of  problems  representing 
various  boundary  conditions  and  end-forces  (including  axial  and  shear  forces) 
because  it  produced  the  largest  discretization  error. 

ELASTIC  RANGE 

Figure  5 indicates  that  the  discretization  error  decreases  monotonical ly 
with  the  element  length.  However,  for  a given  element  length,  the  discreti- 
zation error  increases  with  ec>  This  is  primarily  due  to  the  coupling  of 
axial  and  flexural  deformations,  which  accounts  for  one  source  of  geometric 
nonlinearity  in  the  finite  element  model  (the  formulation  of  conditions  of 
equilibrium  for  the  deformed  state  accounts  for  the  second  source). 

Further  insight  into  the  convergence  characteristics  of  the  finite 

element  model  can  be  gained  by  separating  the  axial  and  flexural  strain 

components  as  indicated  by  Eq.  15.  The  results  are  presented  in  Table  1 for 

three  slenderness  ratios.  They  indicate  that  in  the  limit  e^  ->■  ec  and  sa  > 0, 

which  agrees  with  the  boundary  conditions.  However,  for  a given  element 

length,  the  axial  deformation  term,  £a,  accounts  for  a major  share  of  the 

discretization  error.  For  example,  if  L/r  = 208  and  M = 6 kip  inches,  the 

discretization  error  without  ea  is  2%,  but  with  ea  it  is  14. 7%.  Moreover, the 

c c 

term  ea  grows  rapidly  with  the  flexural  deformations  of  the  finite  element, 

i.e.,  with  the  degree  of  geometric  nonlinearity.  For  example,  if  L/r  = 208 

and  the  moment  increases  from  6 to  36  kip  inches,  ea  increases  relative  to 

e from  16.7%  to  48.9%. 
c 

The  dominant  effect  of  the  term  ca  on  the  strain  state  can  be  explained 
with  the  aid  of  Eq.  2.12  in  reference  1 . For  this  particular  problem,  the 
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Table  1.  - DISCRETIZATION  ERROR 


L/r 

M 

kip  inches 

ccO03) 

rcoo3) 

e^(103) 

ca(103) 

ed 

(%) 

208 

6.0 

-0.150 

-0.128 

-0.153 

0.025 

14.7 

52 

6.0 

-0.150 

-0.148 

-0.150 

0.002 

1 .3 

13 

6.0 

-0.150 

-0.150 

-0.150 

0.000 

0.0 

208 

36.0 

-0.900 

-0.620 

-1 .060 

0.440 

31  .1 

52 

36.0 

-0.900 

-0.843 

-0.903 

0.060 

6.3 

13 

36.0 

-0.900 

-0.896 

-0.900 

0.004 

0.4 

normal  strain  on  the  centroidal  axis  of  the  element  (the  £-axis)  should  be 
zero.  However,  Eq.  2.12  yields: 

e ( ^ ,0 ) = e9  + eC  (18) 

where 

Ea  = (♦ju1  + ^u4)/L  (19) 

and 

eC  = J (^u2/L  + ^u3/L)2  (20) 

It  can  be  seen  from  Eqs.  2.13  - 2 . 1 9 ( i n reference  1)  that  Ea  and  eC  are  poly- 
nomials in  E of  order  1 and  4,  respectively.  Hence,  the  normal  strain  cannot 
vanish  identically  on  the  centroidal  axis  unless  the  two  polynomials  are  of 
the  same  order,  for  instance,  the  order  of  the  polynomial  of  ra  could  be 
increased  from  1 to  4.  But  this  would  require  3 additional  internal  nodes, 
which  would  increase  the  number  of  degrees  of  freedom  of  the  finite  element 
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from  7 to  10.  The  question  of  the  effect  of  additional  internal  nodes  on 
the  efficiency  of  the  finite  element  model  and  the  selection  of  an  optimal 
number  of  internal  nodes  is  under  study.  As  an  alternate  solution,  the 
order  of  the  polynomial  of  eC  could  be  reduced  from  4 to  1 . This  approach 
has  been  adopted  in  the  new  finite  element  model  presented  in  section  V. 


INELASTIC  RANGE 


Figure  6 indicates  that  in  the  inelastic  range  convergence  is  no  longer 
monotone.  Furthermore,  the  discretization  error  does  not  increase  monoto- 
nically  with  the  value  of  the  extreme  fiber  strain  in  compression.  Hence, 
in  the  inelastic  range,  refinement  of  the  finite  element  mesh  no  longer 
guarantees  uniform  improvement  in  response  predictions. 


QUADRATURE  ERROR 


The  Gauss-Legendre  quadrature  is  used  to  evaluate  the  internal  energy 
of  the  finite  element,  which  is  defined  by  the  relation 


1 1 

^ J tt  (c.n)d^dn 


(21) 


where 


-1  -1 


(22) 


★ 

In  Eqs.  21  and  22,  V is  the  volume  of  the  element;  n is  the  internal  energy 
density  at  any  point  U,n)  of  the  longitudinal  plane  of  the  element;  c,n  are 
the  normalized  coordinates  in  the  longitudinal  and  transverse  direction,  res- 
pectively; and  a ,e  denote  stress,  strain  at  any  point  U,n).  Eq.  21  may  be 
resolved  into  two  1 -dimensional  integrals  as  follows: 


> 1 
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and 


= j j "*U)d£ 


(23) 


-1 


1 


(£)  “ ? I tt  (£,n)dn 


i 


-i 


(24) 


where  tt*(?)  is  the  internal  energy  per  unit  length  of  the  element;  A is  the 
cross-sectional  area  and  L is  the  length  of  the  element.  Eqs.  23  and  24 
can  be  expressed  via  the  Gaussian  quadrature  formula  in  the  form 


l n 

* = * Z H.tt*(cJ 
C i=l  1 1 


(25) 


and 


tt*(c)  = * Z H.  ir*(5,r,,) 
£ j=l  J 3 


(26) 


Substitution  of  Eq.  26  into  Eq.  25  yields  the  quadrature  formula  for  Eq.  21 
n m 


w - I I Z H.H  . it*(e.  ,n  .) 

4 i=l  j=l  1 J 1 J 


(27) 


where  H.  and  are  the  weights,  and  £. , n j are  the  Gauss  points  (inter- 
polation points)  in  the  longitudinal  and  transverse  direction,  respectively. 
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EXACT  QUADRATURE 

If  the  internal  energy  density  is  described  by  a polynomial,  one  can 
determine  the  number  of  Gauss  points  required  to  produce  exact  quadratures. 
Suppose  the  constitutive  law  is  expressed  in  polynomial  form  such  that  the 
stress,  a,  is  a polynomial  in  the  strain,  e,  of  order  p.  It  follows  from 
Eq.  22  that  the  internal  energy  density  is  a polynomial  in  e of  order  p+1 . 
Since  the  strain  is  a polynomial  in  £ and  n of  order  4 and  1,  respectively 
(see  Eq.  2.12  in  reference  1 ),  the  strain  energy  density  is  a polynomial 
in  £ and  n of  order  4(p+l)  and  (p+1),  respectively.  Finally,  since  it  takes 
n Gauss  points  to  integrate  a polynomial  of  order  2n-l  exactly,  the  number 
of  Gauss  points  required  for  an  exact  energy  quadrature  is 

n = (4p+5)/2  (28) 

in  the  £-direction  and 

m = (p+2)/2  (29) 

in  the  n-direction,  where  p denotes  the  order  of  the  polynomial  representing 
the  stress  as  a function  of  strain. 

Since  the  strain  field  of  the  new  finite  element  model  (section  V) 
is  described  by  a polynomial  of  order  1 in  £ and  n,  Eq.  29  is  valid  for  both 
directions  of  the  new  finite  element. 

ELASTIC  RANGE 

If  Hooke's  law  is  valid  (i.e.,  if  p=l ) , exact  eneroy  quadratures  can 
be  achieved  with  the  old  finite  element  model  (see  reference  1 ) for  n=5, 
m=2  and  with  the  new  finite  element  model  (see  section  V)  for  n=m=2. 

INELASTIC  RANGE 

Since  the  constitutive  laws  in  SINGER  are  not  described  by  polyno- 
mials (only  by  piecewise  linear  functions),  exact  Gauss  rules  (Gauss  dis- 
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tributions)  cannot  be  defined  in  the  inelastic  range  even  for  monotonic 
loading.  In  any  case,  the  possibility  of  strain  reversals  (hysteretic 
response)  makes  the  polynomial  representation  of  the  complete  constitu- 
tive law  unlikely. 

★ 

Since  the  internal -energy  density  function,  tt  , cannot  be  described 
analytically  in  the  inelastic  range,  formulas  to  estimate  quadrature  errors 
cannot  be  applied  (they  require  knowledge  of  at  least  the  first  order  de- 
rivative of  it  , e.g.,  [18]).  Hence,  error  estimates  can  only  be  obtained 
by  numerical  experimentation.  For  instance,  one  can  establish  an  'optimal' 
Gauss  rule,  defined  as  the  minimum  Gauss  rule  required  to  obtain  the  maxi- 
mum quadrature  accuracy  consistent  with  the  computer  accuracy,  and  assess 
the  quadrature  accuracy  of  coarser  rules  relative  to  the  optimal  rule.  This 
approach  has  been  used  with  the  aid  of  a stream! ined  computer  program  that 
has  variable  Gauss  rules.  (Note  that  presently  the  Gauss  rules  in  SINGER 
are  fixed) . 

LONGITUDINAL  GAUSS  RULE 

Quadrature-discretization  error  studies  of  the  'test'  problem  in  the 
elastic  and  inelastic  range  indicate  that  the  longitudinal  quadrature  error 
corresponding  to  the  3-point  rule  of  SINGER  is  always  significantly  below 
the  discretization  error.  Consequently,  it  appears  that  the  longitudinal 
quadrature  error  can  be  controlled  indirectly  via  the  discretization  error. 
However,  since  it  may  be  more  efficient  to  control  the  two  errors  indepen- 
dently, and  since  this  observation  may  not  always  be  true,  it  is  recom- 
mended that  variable  Gauss  rules  be  introduced  in  SINGER. 

TRANSVERSE  GAUSS  RULE 

A quadrature  error  study  of  the  internal  energy  per  unit  length  as 
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defined  by  Eq.  24  is  based  on  the  reinforced  concrete  section  of  reference  25 
(see  figure  8).  The  purpose  is  to  determine  the  accuracy  of  the  fixed  trans- 
verse rule  in  SINGER  relative  to  the  oDtional  rule  for  a wide  ranqe  of  in- 
elastic deformations  ( see  figures  2 and  6 in  Vol . II  for  a description 
of  the  Gauss  point  locations  used  in  SINGER).  The  results  of  this  study 
indicate  that  the  fixed  rule  in  SINGER  is  quite  accurate  for  predicting 
the  total  energy  stored  in  the  cross-section.  However,  it  may  cause  larne 
errors  in  Dredictinq  the  enerqy  stored  in  the  confined  reqion  of  the  cross- 
section.  For  instance,  for  a qiven  strain  distribution,  the  quadrature  error 
in  the  total  enerqy  is  2.3%,  but  the  error  in  the  enerqy  of  the  confined 
reqion  is  92.8%.  In  spite  of  this  laroe  relative  error,  the  total  Quadra- 
ture error  is  small  since  the  enerqy  contribution  of  the  confined  reaion  to 
the  enerqy  stored  in  the  cross-section  is  only  2.5%. 

Althouqh  the  present  transverse  Gauss  rule  in  SINGER  appears  to  be 
adequate,  the  user  has  no  means  of  estimatinq  or  controllinq  the  resultinq 
errors.  Hence,  it  is  stronqly  recommended  that  a variable  transverse  Gauss 
rule  be  introduced  in  SINGER. 

ERROR  CONTROLS 

It  is  recommended  that  the  initial  selection  of  the  finite  element  mesh 
of  the  structure  and  the  Gauss  rule  of  each  element  be  based  on  a preliminary 
analysis  of  the  "test  problem".  The  cross-sectional  properties  of  the  test 
problem  should  correspond  to  those  of  the  actual  structure  (this  may  require 
several  test  problems  for  one  structure).  The  objective  of  the  preliminary 
analysis  is  to  determine  the  lenqth  and  the  Gauss  rule  of  the  finite  element 
for  each  ranqe  of  response  (e.q.,  elastic  ranqe,  hinqe  element)  that  will 
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keep  the  discretization  and  quadrature  errors  within  specified  bounds. 

The  computation  of  the  exact  value  of  the  extreme  fiber  strain  in 
compression  at  the  fixed  end  can  be  performed  via  classical  methods  or 


directly  with  SINGER.  Since  the  stress  resultant  at  the  fixed  end  is  known, 
the  strain  state  (in  the  elastic  and  inelastic  range)  can  be  found  from  con- 
ditions of  equilibrium.  In  the  inelastic  ranqe  an  iterative  approach  is 
necessary  if  the  strain  state  is  sought  for  a prescribed  moment.  On  the 
other  hand,  one  can  prescribe  the  strain  state  and  compute  the  correspond- 
ing moment  directly  (analogous  to  the  approach  used  to  generate  interaction 
diagrams  for  reinforced  concrete  columns).  Alternatively,  the  strain  state 
at  the  fixed  end  can  be  computed  with  SINGER  to  any  desired  degree  of  ac- 
curacy by  selecting  the  element  length  sufficiently  small. 
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SECTION  V 


NEW  FINITE  ELEMENT  MODEL 

In  order  to  minimize  the  directional  preference  (section  II)  and 
improve  the  accuracy  (section  IV)  of  the  finite  element  model  in  SINGER, 
a new  model  has  been  developed.  It  differs  from  the  old  model  as  follows: 

1.  The  x-axis  of  the  local  moving  frame  reference  is  collinear  with 
chord  that  joins  the  two  external  nodes  of  the  element,  and  the  origin 
of  the  reference  frame  coincides  with  the  midpoint  of  the  chord  (fig. 
7). 

2.  The  strain  field  is  described  by  a polynomial  of  order  1 in  both 
local  coordinates;  i.e.,  the  strain  field  is  linear. 

COORDINATE  TRANSFORMATION 


The  transformation  of  the  global  nodal  displacements  into  local  defor- 
mation components  follows  from  figure  7: 


u1  = a/2 

(30) 

u2  = Ui3  " Y 

(31) 

“3  * Uj3  ■ T 

(32) 

where 

Y = 0 - a 

(33) 

A = | AX* | - | AX | 

(34) 

|AX*|  = [(AX*)2  + (AX*)2]172 

(35) 

|AX|  = [(AX.,)2  + (AX2)2]172  = L 

(36) 

AX  = Xj  - Xi 

(37) 

★ 

AX  = AX  + All 

(38) 
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AU  = IK  - Ui 

(39) 

Y = tan-^ (a/b) 

(40) 

a = AX-jAX*  - AX*AX2 

(41) 

b = aX*aX1  +aX2aX2 

(42) 

In  Eqs.  30  thru  42 

u = one  half  of  the  change  in  element  length 

A = the  change  in  element  length 
u2,u3  = the  relative  element  end-rotations 


U i 3 » U j 3 = rotations  of  the  nodes  i,  j 


IL  ,Uj  = deflection  vectors  of  nodes 


j 


X.,Xj  = position  vectors  of  nodes  i,  j 
Y = chord  rotation 

•ff 

|X  | = chord  length 

L = initial  element  length 

Note  that  the  internal  nodal  deflection,  u^ , represents  a generalized  dis- 
placement of  the  system  model  [1  ];  hence,  no  transformation  is  required. 


FINITE  ELEMENT  MODEL 

The  deformation  field  of  the  discrete  model  is  described  by  four  inter- 
polation functions  corresponding  to  the  four  deformation  components  u-j  , u^ , 
u3,  and  u4  (figure  7).  The  interpolation  functions  are  expressed  in  terms  of 
the  normalized  coordinate 


t = 


-l<t<l 


(43) 


3f 


AXIAL  INTERPOLATION 


I 


On  the  basis  of  the  Lagrange  interpolation  formula  [24],  the  axial 


deflection,  u,  can  be  expressed  as 


u = f i u ( -1 ) + f2u(o)  + ^3uO) 

(44) 

where 

^ = t ( t-1 )/2 

(45) 

f2  = (1+t)  (1-t) 

(46) 

l3  = t(t+l )/2 

(47) 

Imposing  the  continuity  conditions  (figure  7) 

-u(-l)  = u (1 ) = iTj 

(48) 

u(o)  = u4 

(49) 

one  obtains 

u = *1u1  + f4u4 

(50) 

where 

4>1  = t 

*4  = (1-t2) 

(52) 

TRANSVERSE  INTERPOLATION 

With  the  aid  of  the  Hermite  interpolation  formula 

[24],  the  transverse 

deflection,  v,  can  be  expressed  as 

V = f2q2  + f3q3 

(53) 

where 

q2  = u2L/2 

(54) 

q3  = ^3L/2 

(55) 

= (t+1  )(t-l )2/4 

(56) 
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1 


»3  ■ (t-l)(t+l)V4 


Note  that 


clx  $2U2  + *3U3 


where 


= (t-1 ) (3t+l )/4 
^ = (t+1 ) (3t-l )/4 

It  follows  that  v satisfies  the  boundary  conditions  (figure  7) 
v(-l)  = v(l)  = 0 

a?1)- =3 


(57) 

(58) 

(59) 

(60) 

(61) 

(62) 

(63) 


STRAIN  FIELD 

The  strain-displacement  relation  of  the  beam-column  can  be  expressed 
in  the  form  [ 1 ] 


or 


/ \ du  .1  /dv\2  d v 


c\  _ 2 du  . 2 /dvv2  4s  d v 

e(t’s)  - rat + r2(dt}  * r 

L dt 


(64) 


(65) 


where 

s = y/L  (66) 

When  u and  v are  discretized  according  to  Eqs.  50  and  53,  the  first  and  third 
term  on  the  right-hand  side  of  Eg.  65  become  polynomials  in  t of  order  1 
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while  the  second  term  (the  coupling  term)  becomes  a polynomial  in  t of 
order  4.  As  discussed  in  section  IV,  a prerequisite  for  e to  vanish  iden- 
tically along  the  reference  axis  (i.e.,  for  e(t,o)  =0)  is  that  the  poly- 
nomials of  the  first  and  second  term  of  Eq.  65  are  of  the  same  order. 

This  condition  can  be  satisfied  without  increasing  the  degrees  of  freedom 
of  the  element  by  decreasing  the  order  of  the  polynomial  of  the  second 
term  from  4 to  1 . Accordingly,  let 

e<‘-0> ' rar*  p <3r>2  ’ ao + ai  1 <67^ 


where  aQ,  a-|  are  undetermined  coefficients.  The  solution  to  Eq.  67  can  be 


expressed  in  the  form 


f£u(t)-u(-l)>  \ (aQ  + 3lt)dt  - — 2 J (gfrdt 


t 

i ( 

L2  J 


Integration  of  Eq.  68  and  substitution  of  the  continuity  conditions,  Eqs.  48, 
49,  into  Eq.  68  yields  two  algebraic  equations  in  the  unknown  coefficients: 

\ + u4)  = aQ  - i - 2B  (69) 

L U1  ~2ao  " 2A  (70) 

where 


A = 72  i (ar)2  dt 


8 ' M (dar>2  dt 
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The  solution  of  Eqs.  69,  70  yields 

2 - . . 

Vr"]" 

al  ' T “4  + 2(A'2B> 

and  from  Eqs.  67,  73  and  74  one  obtains 

e(t,o)  = - 2tu4)  + (a  + 6t) 

where 

a = A 

0 = 2(A-2B) 


(73) 

(74) 


(75) 


(76) 

(77) 


Note  that  the  first  term  on  the  riqht-hand  side  of  Eq.  75  can  be  obtained 
by  substituting  Eq.  50  into  Eq.  67;  hence,  the  second  term  of  Eq.  75  re- 
presents the  linearized  couplinq  term.  Substitution  of  Eq.  53  into  Eqs. 
71,  72  and  integration  yields  A and  B,  which  can  be  substituted  into  Eqs. 


76,  77  to  yield 

J_  ,-2  1 - - , 2 , 

a '15  u2  ' 2 U2U3  U3  (78) 

» "V6  (=|  - =I>  <79> 

Observe  that  Eq.  65  can  be  expressed  as 

e(t.s)  = c(t,o)  - 

L dr  (80) 

Substitution  of  Eq.  53  into  Eq.  80  yields 

e ( t , S ) = e ( t ,0 ) - s[(3t-1)u2  + (3t+1)u3]  (81) 
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Finally,  on  the  basis  of  Eqs.  75  and  81,  one  can  express  the  strain  field 
of  the  new  finite  element  model  in  the  form 

e(t,s)  = ([-u-j+a)  + (S-fu4)t  - s[(3t-l)u2  + (3t+l)"u3]  (82) 

which  represents  a complete  polynomial  in  t and  s of  order  1 . 
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SECTION  VI 


1 


CAPALi LITY  OF  SINGER 

The  SINGER  program  was  designed  to  predict  the  response  of 
skeletal  structural  systems  incorporating  nonlinear  material  charac- 
teristics as  well  as  finite  displacement  configurations.  The 
system  inertia  forces  included  in  the  solution  process  provide  a 
dynamic  capability.  This  section  is  intended  to  demonstrate  the 
analysis  capability  of  the  current  version  of  the  program  by 
studying  the  type  of  problems  that  can  be  analyzed  as  well  as  the 
accuracy  of  the  macro  response  predictions.  In  addition,  the  con- 
vergence characteristics  of  certain  strain  discontinuities  are 
documented.  In  each  problem,  the  behavior  is  intended  to  repre- 
sent the  entire  continuum  range  up  to  the  ultimate  state.  This 
is  insured  by  supressing  the  local  element  failure  checks  asso- 
ciated with  the  failure  criteria  developed  earlier. 

DEMONSTRATION  PROBLEMS 

The  demonstration  problems  documented  are  primarily  concerned 
with  the  static  response  predictions  of  reinforced  concrete  systems. 
The  two  dynamic  problems  are  single  degree  of  freedom  systems;  one 
problem  has  an  el astopl astic  material,  while  the  other  problem  has 
reinforced  concrete  properties.  Each  problem  is  developed  to 
include  the  description  of  the  system,  materials,  and  the  time 
dependent  forcing  function  if  appropriate.  The  computed  response 
is  shown  in  comparison  with  some  reference  response  data. 

PROBLEM  NO.  1:  REINFORCED  CONCRETE  BEAM 

The  system  analyzed  in  this  problem  is  a simply  supported 
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reinforced  concrete  beam  with  a concentrated  transverse  load 
applied  at  the  center  of  the  span.  The  results  demonstrate  the 
static  prediction  capability  for  a reinforced  concrete  system 
gradually  loaded  to  failure.  The  comparison  basis  is  the  experi- 
mental test  results  of  the  corresponding  physical  beam  measured 
by  the  load-deflection  function  to  failure  published  by  Burns  and 
Siess  [25  ].  Figure  8 shows  the  system  details  and  the  material 
properties;  figure  9 shows  the  load-deflection  response  compari- 
son between  the  computed  and  experimental  values. 

The  comparison  shows  only  small  differences  between  the  two 
functions  up  to  the  limit  point.  The  largest  deviations  occur 
near  the  yield  load  level,  but  limit  point  is  accurately  predicted. 

PROBLEM  NO.  2:  CERF  BEAM  TEST 

This  system  is  a simply  supported  reinforced  concrete  beam 
with  two  equal  lateral  forces  applied  symmetrically  with  respect 
to  the  center  of  the  span.  The  predicted  response  is  intended 
for  comparison  with  specific  experimental  results  obtained  for 
beam  5-0-1  at  the  Civil  Enqineerinq  Research  Facility  [26], 

The  geometry  of  the  beam  is  shown  in  figure  IQ.  The  beam 
was  divided  into  12  elements,  two  of  which  were  steel  end  plates 
which  were  used  for  loading  the  beam  in  the  test  apparatus.  The 
remaining  10  elements  had  cross  sections  as  shown.  The  beam 
element  nodes  were  selected  to  correspond  to  instrumentation  loca- 
tions of  the  experimental  test  beams.  The  material  properties  of 
the  concrete  and  steel  were  chosen  to  match  the  experimental  proD- 
erties;  these  are  shown  in  figure  11  . The  equal  lateral  loads, 
F/2,  were  located  at  nodes  6 and  8. 
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FIG.  10  UNM/CERF  BEAM-COLUMN, 
PROBLEM  NO.  2 


CONCRETE  STRESS-STRAIN 


FIG.  11  MATERIAL  PROPERTIES, 
PROBLEM  NO.  2 
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A comparison  of  the  centerline  transverse  displacements, 
node  7,  from  the  SINGER  results  and  the  CERF  experimental  measure- 
ments is  given  in  figure  12.  The  tensile  steel  strain,  which  is 
the  strain  in  the  bottom  reinforcing  bars,  is  compared  at  node  6 
as  shown  in  figure  13.  The  concrete  strain  at  the  top  of  the  beam 
at  the  centerline,  node  7,  is  also  compared  in  the  same  figure. 

These  results  indicate  that  the  general  behavior  is  accu- 
rately predicted  for  both  the  displacement  measurement  and  the 
strain  measurement  well  beyond  the  linear  response  range.  The 
results  also  show  that  the  computed  response  has  more  stiffness 
in  the  linear  range  than  the  corresponding  experimental  results. 
However,  when  the  displacements  become  large  and  the  beam  is  well 
within  the  inelastic  region  of  response,  SINGER  exhibits  less 
stiffness  than  the  experimental  behavior.  The  SINGER  program 
terminated  execution  when  the  load  was  incremented  beyond  the 
maximum  load  capacity  of  the  beam.  Collapse  conditions  were  not 
accurately  predicted  by  SINGER. 

PROBLEM  NO.  3:  CERF  BEAM-COLUMN  TEST 

The  system  analyzed  for  this  problem  is  identical  to  the 
system  details  shown  in  figure  10.  The  two  equal  lateral  forces 
are  applied  in  combination  with  an  axial  compressive  force.  The 
static  response  prediction  is  intended  for  comparison  with  the 
specific  experimental  results  developed  for  beam  5-2-1  from  CERF  [26l. 

The  material  properties  are  essentially  the  same  as  those 
shown  in  figure  11.  One  difference  is  a lower  crushing  strength 
for  the  concrete  stress-strain  function.  For  this  beam,  the 
crushing  strength  was  approximately  4336.  psi  (29.9  MPa).  The 
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FIG.  12  CENTERLINE  DISPLACEMENT  OF  CERF  BEAM  5-0-1, 
PROBLEM  NO.  2 
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TENSILE  STEEL  STRAIN,  NODE  6 


TOP  FIBER  CONCRETE  STRAIN,  NODE  7 


FIG.  13  BEAM  5-0-1  STRAIN  COMPARISONS,  PROBLEM  NO.  2 
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loading  consisted  of  the  two  lateral  forces  F/2  ar.d  the  combined 
axial  force  P.  The  ratio  P/F  of  1.91  was  developed  for  each  load 
increment  to  failure. 

The  response  comparison  is  based  on  displacement  and  strain 
measurements.  The  centerline  displacement,  node  7,  for  both  SINGER 
and  the  experimental  test  measurements  is  given  in  figure  14. 

A comparison  of  the  tensile  strains  in  the  bottom  reinforcing  steel 
at  node  6 is  shown  in  figure  15.  The  concrete  strain  at  the  top  of 
the  beam  at  the  centerline  (node  7)  is  also  illustrated  in  figure  15. 
The  computed  results  accurately  predict  the  initial  response  and 
the  limit  point  but  Fail  to  define  the  gradual  loss  of  stiffness 
near  the  yield  load.  Tne  decreasing  experimental  load-deflection 
path  beyond  the  limit  point  cannot  be  developed  by  the  program 
since  the  equilibrium  path  beyond  the  limit  state  is  unstable. 
Experimentally,  this  region  of  the  function  can  be  defined  since 
the  displacements  can  be  physically  controlled  and  the  forces 
adjusted  accordingly.  Comparison  of  centerline  displacement  quan- 
tities of  beam  5-2-1  indicates  that  the  finite  element  model  of 
SINGER  is  stiffer  than  the  experimental  beam  all  the  way  to  col- 
lapse. The  tensile  steel  strain  comparison  reflects  this  behavior. 
This  demonstrates  that  the  steel  response  may  dominate  the  system 
response. 

PROBLEM  NO.  4:  REINFORCED  CONCRETE  WITH  CYCLIC  AXIAL  LOAD 

The  system  analyzed  is  a reinforced  concrete  member  with  a 
force  applied  parallel  to  the  longitudinal  axis.  The  force  is 
first  incremented  as  a compression  force  up  to  the  value  of 
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FIG.  14  CENTERLINE  DISPLACEMENT  OF  C ERF  BEAM  5-2-1, 
PROBLEM  NO.  3 
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FIG.  15  BEAM  5-2-1  STRAIN  COMPARISONS,  PROBLEM  NO.  3 
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650. X10  lbs.  so  that  hath  the  steel  and  concrete  are  active.  From 

this  point,  the  force  is  reduced  to  zero,  then  applied  as  a tension 
3 

force  up  to  150. X10  lbs.  At  this  value  it  is  again  reduced  to 

3 

zero  and  then  applied  as  a compression  force  up  to  150. X10  lbs. 

The  magnitudes  reached  by  the  applied  force  are  sufficiently  large 
to  produce  inelastic  response,  including  strain  hardening  in  the 
steel  during  the  tension  range.  The  purpose  of  the  problem  is  to 
demonstrate  the  accuracy  of  the  static  inelastic  response  of  a 
reinforced  concrete  member  subjected  to  both  compression  and  tension. 
Since  the  resulting  strain  states  are  uniform,  the  response  can  be 
checked  by  exact  calculations  for  the  axial  displacement  and  the 
corresponding  axial  force  at  an  equilibrium  configuration.  The 
system  and  materials  used  are  shown  in  fiqure  16;  the  response  and 
comparison  checks  are  shown  in  figure  17.  The  check  values  agree 
perfectly  with  the  computed  response.  Very  good  convergence  charac- 
teristics are  usually  obtained  for  the  solution  of  axially  loaded 
systems  for  the  complete  range  of  response.  The  large  displacement 
values  in  the  tension  range  are  due  to  the  lack  of  concrete  strength 
in  tension.  The  tensile  force  is  resisted  only  by  the  reinforcing 
bars. 

PROBLEM  NO.  5:  SINGLE  DEGREE  OF  FREEDOM  ELASTIC-PLASTIC  SPRING 

This  system  consists  of  an  ideal  massless  spring  (bar)  with 
mass  attached  at  the  free  end.  A time-dependent  forcing  function  is 
applied  to  the  mass  in  the  axial  direction  so  that  the  resulting 
system  has  only  a single  degree  of  freedom.  The  material  has 
ideally  elastic-plastic  properties  with  a linear  unloading 
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charateristi c;  the  forcing  function  is  a decreasing  triangular 
pulse.  The  system,  material  function,  and  forcing  function  are 
shown  in  figure  18. 

The  purpose  is  to  demonstrate  the  dynamic  capability  with 
inelastic  response  using  a problem  which  can  be  checked  by  exact 
computations.  The  displacement  response  and  comparison  checks 
are  shown  in  figure  19.  Although  the  computed  and  exact  functions 
are  almost  identical  in  the  initial  portion  of  the  response,  accu- 
mulated errors  eventually  cause  some  discrepancy  to  appear  and 
grow.  The  error  observed  in  this  problem  is  primarily  a phase 
shift  in  the  displacement  function. 

PROBLEM  NO.  6:  REINFORCED  CONCRETE  BEAM  WITH  AXIAL  FORCE  PULSE 

This  problem  analysis  develops  the  response  for  a massless 
reinforced  concrete  element  with  attached  mass  at  the  free  end.  The 
single  degree  of  freedom  system  has  an  applied  force  pulse  in  the 
axial  direction.  The  material  properties  are  linearized  so  that 
a closed  form  solution  can  be  used  as  a comparison  basis.  The 
element  details  and  the  forcing  function  are  shown  in  figure  20. 

Both  concrete  and  steel  properties  are  defined  in  terms  of 
their  respective  modulus  of  elasticity.  For  concrete,  the  elastic 
modulus  is  34.474  GPa  , while  the  steel  modulus  is  241.315  GPa 
The  forcing  function  is  a linearly  decreasing  triangular  force 
pulse  which  varies  from  650  kN  at  time  zero  to  a zero  value  at 
time  0.002  second  . 

The  closed  form  displacement  response  during  the  loading  phase  is 
U(t)  = P0/k(sin(u)t)/(1)t1  - cos(ut)  - t/tj  + 1) 
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where  k is  the  stiffness  defined  as  AE/L;  A is  the  transformed  cross- 
sectional  area  of  concrete,  E is  the  elastic  modulus  of  concrete  and 
L is  the  length  of  the  element.  The  natural  frequency,  w is  defined 
as  /(k/M  , where  M is  the  concentrated  x-component  of  mass.  The 
constant  t^  is  the  time  at  which  the  force  has  decreased  to  zero, 
0.002  second. 

Comparison  of  the  SINGER  response  with  the  theoretical  response 
is  shown  in  figures  21  and  22.  The  comparison  is  only  made  through 
the  point  where  the  beam  is  still  in  compression.  The  closed  form 
solution  would  have  to  be  re-evaluated  for  the  response  when  tension 
occurs  since  the  stiffness  will  be  changed.  For  the  time  of  response 
that  is  compared,  there  is  negligible  difference  between  the  SINGER 
prediction  and  the  closed  form  solution.  When  the  beam  goes  into 
tension,  the  response  time  increases  indicative  of  the  longer  period. 
CONVERGENCE  CHARACTERISTICS  OF  STRAIN  DISCONTINUITY 

It  has  been  observed  in  computed  output  that  a single  element 
required  to  represent  relatively  large  changes  in  strain  within 
its  length  may  in  fact  create  a strain  discontinuity  with  respect 
to  an  adjacent  element  while  achieving  a system  equilibrium  configu- 
ration. Large  changes  in  strain  within  an  element  occur  when 
one  element  is  severely  strained  beyond  the  yield  state  of  the 
material  while  an  adjacent  element  is  much  less  severely  yielded, 
or  it  is  not  yielded  at  all.  The  strain  values  at  the  sections  on 
either  side  of  a connecting  joint  may  differ  not  only  in  magnitude 
but  also  in  sign.  This  discontinuous  behavior  has  been  observed 
in  computed  results  for  reinforced  concrete  elements  as  well 
as  for  steel  wide  flange  elements.  The  number  of  elements 


63 


64 


FIG.  22  COMPARISON  OF  VELOCITY  AND  ACCELERATION 
OF  CONCRETE  BEAM,  PROBLEM  NO.  6 
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chosen  to  represent  the  system  can  significantly  affect  the  form 
and  severity  of  the  discontinuities  created.  The  strain  values  which 
occur  as  a part  of  the  output  record  for  measuring  the  interelement 
strain  continuity  are  those  values  computed  at  the  points  of  the  end 
sections  of  each  element  after  equilibrium  has  been  achieved. 

The  characteristics  of  strain  discontinuity  are  documented  in 

figures  23  and  24  with  the  specific  objective  to  show  the  relative 

strain  magnitudes  which  occur  and  to  demonstrate  the  effect  of 

increasing  the  number  of  elements.  The  problem  data  are  identical  to 

the  data  for  problem  No.  1 (figure  8)  with  slight  differences  in 

the  material  properties:  f'  = 4000.  psi  anJ  f = 40,000.  psi . The 

c y 

strain  in  the  tension  steel  is  used  to  measure  the  strain  behavior. 
Figure  23  shows  the  strain  distribution  at  the  first  load  level  to 
cause  the  discontinuous  behavior;  figure  24  shows  the  distribution 
at  the  next  load  increment. 

In  a reinforced  concrete  section,  the  tension  requirement  for 
bending  is  satisfied  by  the  reinforcing  steel  only;  in  SINGER,  the 
concrete  material  has  no  tensile  strength.  The  effect  of  this  prop- 
erty is  that  once  the  steel  begins  to  yield,  the  strain  passes 
almost  immediately  into  the  strain  hardening  range  since  there  is 
no  other  source  of  positive  stiffness  to  resist  increased  loads. 

The  strain  is,  therefore,  either  below  the  yield  strain  or  into 
the  strain  hardening  region.  This  severe  strain  qradient  near  the 
section  of  maximum  bending  moment  corresponds  to  larqe  deformation 
changes  in  a highly  localized  region  of  the  beam.  As  shown  in 
figure  23,  the  strains  change  from  yield  strain  (0.00138)  to  the 
order  of  twice  strain  hardening  (0.023)  within  a distance  of  0.9 
times  the  depth  of  the  beam.  In  figure  24,  the  next  higher  load 
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level,  the  strains  change  from  yield  to  three  times  strain  hardening 
within  a distance  of  1.35  times  the  depth  of  the  beam.  The  average 
rate  of  change  for  both  regions  is  essentially  the  same.  The  strains 
within  the  remaining  portion  of  the  beam  are  at  or  below  yield  strain. 

These  strain  gradients  must  be  modeled  by  the  SINGER  elements 
which  are  constrained  to  produce  linear  strain  fields.  It  has  been 
previously  discussed  (see  section  IV)  that  the  continuum  strain 
distribution  in  an  element  becomes  complex  and  discretization  errors 
grow  when  material  properties  are  nonlinear  and  deformations  are 
large.  The  lack  of  modeling  capability  which  would  provide  for 
continuous  strain  distribution  in  the  elements  is  reflected  in  the 
discontinuous  characteristics  that  occur  in  these  complex  strain 
regions.  The  strain  continuity  can  be  maintained  by  elements  if 
the  continuum  strain  field  can  be  closely  approximated  by  a linear 
distribution.  This  is  demonstrated  in  the  less  severely  strained 
regions  of  the  beam  represented  in  figures  23  and  24  . 

The  computational  reason  for  the  possibility  of  the  discon- 
tinuity in  strains  at  adjacent  sections  is  that  the  minimization 
process  does  not  include  any  direct  control  on  strain  values.  The 
generalized  coordinate  values  corresponding  to  an  equilibrium 
configuration  are  determined  by  the  minimization  of  the  total 
system  energy.  The  total  energy  is  composed  of  the  summation  of 
the  individual  element  energies  and  the  contribution  from  the 
applied  forces  and  the  inertia  forces.  The  convergence  of  the 
minimization  process  is  controlled  by  a measure  of  the  largest 
stepsize  used  in  the  search  at  a point  in  the  configuration  space. 

The  error  associated  with  a converged  state  is  determined  by  a 
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dimensionless  measure  of  the  largest  value  of  the  residual  forces 
at  the  joints.  Both  of  these  process  measurements  are  macro 
quantities.  There  is  no  direct  measurement  of  strain  used  to 
control  the  process  or  to  measure  errors.  Its  effect  enters  only 
indirectly  through  the  element  energy  computations.  Hence,  strain 
continuity  between  elements  is  possible,  but  certainly  not 
guaranteed  by  the  process  itself. 

SHEAR  EFFECTS  IN  REINFORCED  CONCRETE 

The  complete  planar  response  of  a beam-column  member  includes  axial, 
flexural  and  shear  deformation  components,  with  the  flexural  and  shear 
deformation  comorisino  the  transverse  displacement.  For  most  commonly 
used  dimensions,  the  flexural  contribution  dominates  the  tranverse  dis- 
placement. However,  it  is  commonly  known  from  the  analysis  of  two  di- 
mensional stress  fields  that  members  with  relatively  laroe  depth  - lenqth 
ratios  (deep  beams)  may  have  a significant  part  of  the  transverse  dis- 
placement contributed  by  the  shear  deformation  component  of  the  response. 

In  addition,  the  total  shear  deformation  in  reinforced  concrete  deep 
beam  members  loaded  to  failure  may  be  affected  by  internal  concrete 
destruction  and  slippaqe  associated  with  the  qrowth  of  inclined  cracks, 
since  the  failure  state  may  be  significantly  laroer  than  the  initial 
crackino  load  [27]. 

The  element  model  in  SINGER  includes  axial  and  flexural  deformation 
components,  but  it  does  not  include  shear  deformation.  C.onseouentl v , the 
response  of  any  member  is  predicted  accordinq  to  the  deformation  coordinates 
u^  , u,,,  u^,  and  u^  (see  fiqure  7).  For  the  line  element  used  in  SINGER, 
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shear  deformation  could  be  included  by  an  indirect  procedure;  i.e.,  the 
shear  deformation  could  be  annroximated  from  an  averane  shear  strain 
distribution  at  a section  corresDOndino  to  a shear  force  which  is  re- 
quired for  equilibrium  of  the  forces  actinq  on  the  element.  This  ap- 
proximation could  also  be  incorporated  in  the  transverse  deformation 
function.  However,  the  resultinq  expression  is  based  on  the  linear 
material  ranae.  The  extension  into  the  inelastic  material  response  plus 
the  absence  of  test  data  for  this  effect  in  reinforced  concrete  create 
sufficient  uncertainties  so  that  the  resultinq  shear  deformation  pre- 
diction was  judqed  to  be  unreliable  compared  to  the  other  deformation 
components  in  the  element. 

The  effect  of  shear  on  the  failure  of  an  element  is  included  in  the 
failure  criteria  by  an  indirect  procedure  [1].  The  purpose  is  to  de- 
termine the  force  conditions  which  most  probably  cause  a diaqonal  crack- 
inq  failure  state  within  an  element.  It  is  an  indirect  procedure  based 
on  experimental  data,  since  there  is  no  way  to  relate  the  shear  force  at 
a section  with  the  correspondinq  strain  state  of  the  element.  These  criteria 
also  include  the  effect  of  web  reinforcement. 

To  demonstrate  that  the  SINGER  element  model  is  not  capable  of  ac- 
curately predictino  the  response  of  members  which  have  a siqnificant  shear 
deformation,  four  deep  beam  test  problems  were  studied.  The  comparison 
basis  for  the  predicted  response  was  obtained  from  experimental  test 
results  published  by  dePaiva  and  Siess  \27~\.  The  four  problems  correspond 
to  their  beam  identifications  noted  G33S-11 , G33S-12.  G33S-31 , and  G33S-32. 
The  beams  G33S-11,  12  are  referred  to  as  beam  (a.),  and  G33S-31 , 32  are 
referred  to  as  beam  (b.).  The  physical  details  are  summarized  in  finure 
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25;  the  response  comparisons  are  clotted  in  figures  26,  27  and  28.  The 
functional  comparison  relates  the  total  aoplied  load  to  the  central  trans- 
verse deflection  for  each  beam  to  failure,  shown  in  fiaures  26  and  27. 
Figure  28  shows  a comparison  of  the  strain  at  the  center  of  the  tensile 
reinforcing  steel.  In  all  cases  it  is  clear  that  the  SINGER  element  is 
much  stiffer  (less  deformation  at  the  same  load  level).  Since  the  SINGER 
element  includes  only  the  flexural  transverse  deflection,  a significant 
shear  deformation  is  indicated  in  the  test  results  (more  than  50  percent 
of  the  total  transverse  deflection  for  all  four  cases).  The  central  strain 
values  of  fiaure  28  support  these  observations,  but  the  difference  between 
the  predicted  and  exoerimental  values  are  relatively  smaller,  indicating 
closer  agreement. 

The  implementation  of  the  SINGER  failure  criteria  for  these  problems 
Droduced  predictions  of  the  shear  effects  on  element  failure.  These  re- 
sults were  comoared  with  the  exoerimental  inclined  crackino  load,  shown 
below. 


Beam 

Wi  th 

Load  at  Inclined 
Web  Reinf. 

Cracking  (Kips) 

Without  Web  Reinf. 

SINGER 

Experimental 

SINGER 

Experimental 

(a.) 

20.0 

12.90 

2.0 

13.32 

(b.) 

13.0 

16.96 

1 .0 

16.96 

In  SINGER,  the  initial  detection  of  an  inclined  crack  in  a element  without 
web  reinforcement  is  indicated  as  a shear  failure.  For  an  element  with 
web  reinforcement,  the  failure  is  indicated  after  the  web  reinforcement 
yields.  The  experimental  values  tabulated  above  do  not  indicate  failure; 
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Beam 

Material 

Property 

Stress  Value  (ksi) 

With  Web  Reinf. 

Without  Web  Reinf. 

(a.) 

f 1 ! 

c 

2.89 

3.38 

fy  (tens. ) 

47.3 

47.3 

f (compr. ) ; 

51.5 

51.5 

(b.) 

V 

2.91 

2.89 

f (tens. ) 

44.2 

45.2 

f (compr.) 

50.3 

50.3 

f 1 = concrete  strength  in  compression 

f (tens.)  = yield  stress  of  tension  steel 
f (compr. ) = yield  stress  of  compression  steel 


FIG.  25  SYSTEM  AND  MATERIAL  PROPERTIES 
FOR  SHEAR  EFFECTS 
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Computed 


Experimental 


With  Web  Reinforcement 
Without  Wet  Reinforcement 


1 


FIG.  28  TENSION  STEEL  STRAIN  COMPARISON  - BEAM  (b.) 
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they  do  indicate  the  load  at  the  detection  of  the  first  inclined  crack. 
For  the  case  of  no  web  reinforcement,  aqreement  is  poor.  The  aqreement 
for  the  case  with  web  reinforcement  is  reasonable. 

It  should  be  noted  that  a thorouqh  study  of  the  failure  criteria  has 
not  been  attempted  in  this  report.  Such  a study  would  be  an  extensive 
undertakinq  in  itself,  with  a detailed  knowledqe  of  the  capability  and 
accuracy  of  the  basic  proqram  serving  as  a mandatory  prerequisite.  The 
report  has  furnished  the  details  of  the  refinement  and  demonstration  of 
the  current  capability  of  SINGER. 
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PROGRAM  MODIFICATIONS  AND  SUGGESTED  IMPROVEMENTS 

The  objective  of  this  section  is  to  document  the  significant 
modifications  incorporated  in  the  original  version  of  the  SINGER 
program  [28]  and  to  discuss  improvements  which  should  be  made  to 
provide  more  accurate  and  efficient  response  predictions. 

PROGRAM  MODIFICATIONS 

The  original  version  of  the  SINGER  program  was  unable  to  solve 
a significant  class  of  structural  problems  consistent  with  its  potential 
capability,  particularly  those  with  reinforced  concrete  elements. 

This  deficiency  was  due  to  several  factors: 

a.  The  scope  of  the  program  was  extremely  large, 
requiring  many  computational  checks  and  much  time 
to  verify  that  all  parts  were  working  together 
correctly; 

b.  The  time  constraints  left  some  of  the  subroutines 
virtually  unchecked  as  computational  units  (primarily 
those  involved  with  the  internal  energy  computations 
for  reinforced  concrete  elements); 

c.  A complex  storage  scheme  of  selected  arrays,  althouqh 

| 

efficient  for  computations,  was  difficult  to  implement 
and  check; 

d.  Many  people  with  varying  programming  experience 
contributed  to  the  original  program,  creating  a source 
of  coding  errors  and  inefficiency. 
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A considerable  effort  was  required  to  expand  the  usable  capacity 
of  this  program.  The  development  of  this  task  was  influenced  by  the 
evolutionary  growth  of  experience  with  the  code.  The  first  task 
involved  removing  coding  errors  and  undefined  variables  so  that  the 
program  would  execute  a set  of  input  data.  Once  this  was  achieved, 
some  modifications  could  be  made  to  improve  its  efficiency  and  accuracy. 
There  was  no  clear  demarcation  between  these  tasks,  and  an  interaction 
was  unavoidable. 

The  modifications  made  to  the  program  cover  a broad  scope.  The 
most  significant  contributions  are  listed  in  three  categories:  (1) 
correction  of  coding  errors  and  computational  errors;(2)  improvements 
in  computational  efficiency;  and(3)  improvements  in  solution  accuracy. 
The  name  in  parentheses  refers  to  the  corresponding  subroutine  where 
the  modification  was  made. 

1.  Correction  of  coding  errors  and  computational  errors: 

a.  At  least  one  undefined  variable  was  detected  in  each 

of  10  subroutines.  In  many  cases  these  were  corrected 
by  adding  COMMON  statements.  (ADYN,  BEAM,  BOND,  FAIL, 

FORK,  GIDE,  LEAF,  OUTS,  SECT,  STEN). 

b.  The  subtraction  of  steel  reinforcing  bar  area  from 
concrete  area  was  done  incorrectly  (COEN). 

c.  The  accumulation  of  the  energy  contribution  from  the 
reinforcing  bars  was  done  incorrectly  (STEN). 

d.  The  default  functions  defining  the  stress-strain 
properties  for  both  confined  and  unconfined  concrete 
had  some  of  the  constant  values  incorrectly  defined. 
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or  they  were  inconsistent  with  the  User’s  Guide  (MATP). 

e.  A subscript  computation  for  reference  to  a value  stored 
in  DATA  was  computed  incorrectly.  This  particular  error 
caused  serious  computational  problems  (COEN). 

f.  The  computation  for  tne  energy  associated  with  the  steel 
reinforcing  bars  was  incorrectly  done  due  to  an  error 

in  the  Gaussian  quadrature  equations  (STEN). 

g.  A subscript  for  the  DATA  array  was  incorrectly  computed 
for  use  in  generating  forcing  function  values.  An  error 
in  the  forcing  function  value  was  caused  for  a forcing 
function  defined  by  only  two  points  (TABL). 

h.  An  error  occurred  in  the  energy  density  computation 
for  concrete  for  the  special  case  when  the  previous 
and  current  strain  values  were  identical.  This  error 
was  also  of  a serious  nature  preventing  correct 
convergence  in  the  minimization  pru^ess  (CRET). 

2.  Improvements  in  computational  efficiency: 

a.  Several  logical  IF  decisions  were  rewritten  to  avoid 
awkward  GO  TO  transfers. 

b.  Material  data  co  putations  were  removed  from  a subroutine 
within  the  solution  process  to  avoid  unnecessary  repetition 
in  the  computations  which  were  actually  required  to  be 
done  only  oncp  at  the  time  of  data  input  (COEN). 

c.  Many  of  the  subroutines  contained  unreferenced  statement 
numbers  and  redundant  or  duplicate  statements,  which  were 
corrected. 
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d.  Some  segments  were  rewritten  to  use  fewer  statements 
to  accomplish  the  same  objective,  particularly  those 
using  logical  IF  statements  in  a decision  sequence. 

e.  At  least  one  major  subroutine  was  rewritten  to  avoid 
unnecessarily  complicated  organization,  which  was 
difficult  to  follow  (ASAN). 

f.  Revisions  were  made  in  a few  cases  so  that  subroutines 
were  called  only  when  necessary  for  the  given  conditions 
rather  than  having  the  subroutine  make  the  decision 
after  being  called. 

3.  Improvements  in  solution  accuracy: 

Two  major  revisions  were  implemented  with  the  objective 
to  remove  or  reduce  the  directional  properties  of  the  original 
element  model. 

a.  The  first  revision  reduced  the  axial  strain  contribution 
of  the  coupling  term  in  the  strain-displacement  relation 
from  a fourth  order  function  in  the  x coordinate  to  a 
linear  form  to  be  consistent  with  the  linear  flexural 
strain  contribution. 

b.  The  second  revision  (documented  in  section  V)  completely 
redefined  the  basis  for  the  element  reference  coordinate 
system.  The  deformation  coordinates  were  redefined  and 
the  axial  component  of  strain  was  again  linearized  with 
respect  to  the  x coordinate. 

Although  both  modifications  were  useful  in  improving  the  performance 
of  the  element  within  the  system,  the  latter  became  a permanent 
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contribution  to  the  program.  The  improvement  in  the  element  behavior 
was  measured  by  observing  the  symmetry  of  element  deformation  and 
strains  for  a problem  whose  exact  response  had  symmetric  properties. 

SUGGESTED  IMPROVEMENTS 

The  current  version  of  the  SINGER  code,  although  not  yet  completely 
operational,  is  capable  of  solving  a broader  class  of  problems  than 
the  original  version,  including  reinforced  concrete  systems.  The 
progress  made  in  understanding  and  improving  the  program  has  also 
generated  many  ideas  for  further  improvement.  While  some  of  these 
improvements  are  simply  revisions  to  the  form  of  output  for  a problem, 
others  involve  profound  changes  which  would  require  extensive  study 

before  implementation.  The  ideas  considered  to  be  the  most  important 

¥ 

are  discussed  below  under  four  headings:  (1)  detection  of  errors  and 
computational  difficulties-^  2)  improvements  in  the  form  of  input  and 
output;(3)  improvements  in  computational  efficiency;  and(4)  improve- 
ments in  solution  accuracy. 

1.  Detection  of  errors  and  computational  difficulties: 

a.  The  program  should  be  capable  of  computing  and 
displaying  the  components  of  the  element  and 
system  energy  quantities  corresponding  to  a 
prescribed  displacement  configuration  without 
including  the  minimization  process  for  equilibrium. 

This  would  be  useful  for  making  checks  on  the 
internal  energy  computations. 

b.  The  default  unconfined  concrete  material  function 
should  be  revised  to  include  a large  strain  value 
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corresponding  to  a small  stress  value  at  material 
point  7,  similar  to  the  form  used  for  confined  concrete. 

The  area  under  the  function  between  points  6 and  7 could 
be  made  sufficiently  small  to  eliminate  a significant 
energy  contribution.  The  purpose  would  be  to  avoid 
the  possibility  of  strain  values  exceeding  the  final 
material  point  during  the  solution  process. 

c.  There  should  also  be  a check  on  user  input  material 
functions  to  prevent  strains  from  exceeding  the  maximum 
strain  point  during  a problem  solution. 

d.  It  may  be  useful  to  provide  an  option  for  implementing 
the  failure  criteria  for  a specific  problem.  Some  forms 
of  behavior  may  be  more  effectively  studied  by  tracing 
the  system  response  without  incorporating  local  failure 
termination,  and  it  may  be  helpful  for  checking  the 
details  of  the  failure  criteria.  Other  nonrelated 

but  essential  computations  performed  in  the  subroutine 
(FAIL),  such  as  stress  resultant  computations,  should 
remain  operational. 

e.  The  subroutine  which  computes  the  internal  energy  in 
leaf  spring  elements  (LEAF)  is  not  consistent  with  the 
latest  element  model  because  the  element  deformation 
components  have  been  redefined.  This  subroutine,  if 
it  is  to  be  used,  should  be  rewritten. 

f.  It  may  be  important  to  have  a provision  to  dump  selected 
regions  of  the  DATA  and  KDATA  arrays  to  be  able  to  study 
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their  contents  in  conjunction  with  certain  computational 
errors  detected.  This  storage  scheme  has  been  the  source 
of  past  errors  and  it  is  difficult  to  monitor. 

2.  Improvements  in  the  form  of  input  and  output: 

a.  The  default  confined  concrete  stress-strain  values 
which  are  affected  by  the  confinement  factor  Uc(5), 
f c ( 5 ) and  cc(6))  are  never  printed  as  a part  of  the 
output  record.  Their  values  are  not  known  at  the  time 
of  the  function  generation  since  the  confinement  factor 
depends  on  stirrup  spacing,  which  is  input  after  the 
material  properties.  Some  provision  should  be  made  to 
print  the  complete  function. 

b.  Related  to  a.  above,  the  user  input  function  for  confined 
concrete  should  be  printed  out  in  its  entirety  since 

it  is  not  affected  by  the  confinement  property.  Presently, 
the  same  three  values  are  masked  in  the  output  as  it  is 
done  for  the  default  function  since  both  functions  are 
printed  by  the  same  statement. 

c.  The  complete  applied  force  vector  should  be  a part  of  the 
output  record  for  each  load  increment  regardless  of  the 
print  level  desired.  In  the  present  form,  these  forces 
are  not  printed  for  certain  print  level  options. 

d.  It  would  be  desirable  to  include  the  option  to  print 
out  results  in  a dynamic  analysis  only  at  specified 
time  points  to  avoid  extremely  large  volumes  of  output. 
Short  messages  could  be  printed  between  major  output 
points  to  monitor  the  significant  features  of  behavior. 
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3.  Improvements  in  computational  efficiency: 

a.  The  confinement  factor  is  computed  for  every  case  which 
includes  lateral  reinforcement.  However,  it  is  used 
only  for  the  default  confined  concrete  function.  It 
does  not  need  to  be  computed  for  user  input  functions 
because  their  values  are  not  affected  by  this  factor. 

b.  It  could  be  helpful  to  develop  a way  to  monitor  the 
system  response  to  provide  information  related  to  how 
close  the  system  is  to  the  conditions  of  the  limit 
point. 

4.  Improvements  in  solution  accuracy: 

a.  A variable  Gauss  point  mesh  in  an  element  would  be  an 
important  feature  to  control  the  accuracy  in  the  element 
energy  computations  (see  section  IV). 

b.  It  is  necessary  to  study  the  behavior  of  higher  order 
elements  to  be  able  to  decide  on  the  most  efficient 
model  for  computing  the  system  energy  in  the  range  of 
complex  behavior,  with  particular  emphasis  on  improving 
the  interelement  strain  discontinuity  problem  (see 
section  VI ) . 

c.  It  would  be  desirable  to  provide  a way  to  monitor  the 
change  in  slope  for  the  individual  elements  and  perhaps 
impose  a slope  change  limitation  based  on  the  accuracy 
of  element  strain  computations.  It  would  provide  the 
user  with  a warning  when  individual  elements  were 
excessively  deformed.  This  monitor  could  also  include 
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a check  on  the  element  strain  gradients,  and  ultimately 
it  could  be  related  to  a warning  of  possible  strain 
discontinuities. 


The  experience  gained  by  working  with  the  SINGER  program  has 
established  a firm  confidence  in  its  ability  to  solve  problems  within 
its  intended  capability.  Although  this  potential  has  not  been  fully 
realized  to  date,  it  can  ultimately  be  developed. 
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ERRATA  TD  ORIGINAL  REPORT 

The  work  function,  W,  defined  by  Eq.  3.14  in  reference  1 should 
be  replaced  by  the  energy  function,  E,  which  is  defined  as 


n r6m.  x,  . 

E = l ~T  ^ ' aixbi } - (fKi  + Oxh 

i=lLAt^  ^ 1 


bi  bi ' bi 
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+ U 


0) 


where 


At 

a-  = X ■ + X .At  + X,.  —5- 

i ai  ai  ai  3 


The  energy  function  assumes  a stationary  value  if 


(2) 


3E 


8Xbi 


= 0,  i = 1 ,2, . . . ,n 


Eqs.  3 yield  n equations  of  equilibrium  of  the  form 

Vbi  + % = fbi  + f6bi’  1 = 


(3) 


(4) 


where 


xbi  .2  ^xbi  " “i ^ 

At 

Eq.  5 follows  from  Eqs.  3.4  and  3.5  in  reference  1 with  t = At;  i.e.. 


(5) 


xbi  = xi(At^ 
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